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ABSTRACT 

The  methods  of  Quasilinearization,  Dynamic  Programming,  and 
Invariant  Imbedding  are  of  great  practical  use  in  transforming  two 
point  boundary  value  problems  into  forms  that  are  more  readily  solved 
numerically,  i'lach  of  the  three  methods  is  discussed  through  the  use  of 
an  example. 


THREE  HEW  COMPUTATIONAL  METHODS  FOR  SOLVING 


TWO  POINT  BOUNDARY  VALUE  PROBLEMS 

* 

Kenneth  B.  Bley 
The  RAND  Corporation 


The  desirable  goal  of  optimization  often  leads  to  the  undesirable 
requirement.  for  the  solution  of  two  point  boundary  value  problems.  This 
paper  dis  usses,  through  an  example,  three  i, methods  which  reduce  the  two 
point  boundary  value  problem  to  a  more  tractable  form  suitable  for 
solution  on  digital  computers 

To  illustrate  these  methods  we  will  consider  a  single  cost  fun  - 
tion  which  we  will  seek  to  optimize  (i.e.  minimize)  using  each  of  the 
three  techniques. 

Suppose  we  wish  to  find  the  motion  of  a  particle  in  an  inverse 
square  gravity  field.  Hamilton ' s  principle  states  that  the  motion  is 
given  by  the  function  that  minimizes  the  time  integra-  of  the  Lagrangian 
where  the  initial  and  final  states  of  the  particle  are  given.  If  an 
initial  state  and  a  time  are  specified  then  we  have  a  free  boundary  prob¬ 
lem  and  the  velocity  of  the  particle  must  be  zero  at  the  specified  time 
We  therefore  define  a  cost  functional  j[u]  as 


J[u]  3 1  (i  *  *  i)dt 

where  we  have,  as  boundary  conditions, 

u(o)  ■  c 
u(T)  3  o 


(1) 


(2) 


* 
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1.  EULER  EQUATIONS  AND  QUASILINEARIZATION 

In  the  proceeding  section  we  defined  a  cost  function  which  we 
wish  to  minimize  through  the  choice  of  a  suitable  u(t).  If  we  write 

T 

J  u]  ■  F(u,u)dt  (3) 

b 


then  Euler's  Equation, 


^u  dt 


(V\ 


(«0 


which  is  a  necessary  condition  for  an  extremum  of  the  integral,  leads 
to  the  nonlinear  second  order  differential  equation 


...  A. 
u  +  -j  -  0 

u 


(5) 


We  will  reduce  Eq  (5)  to  two  first  order  equations  by  taking 


u  ■  v 


g(u,v) 


V  ■  -  ~  -  h(u, v) 
u 


(6) 


subject  to  the  boundary  conditions 

u(o)  ■  c 
v(T)  ■  o 


(7) 


We  are  now  ready  to  use  the  method  of  Quasilinearization  to  solve 
Eq  (6).  To  apply  this  method  we  will  expand  Eq  (6)  in  a  Taylor's 
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series  in  the  functions  u  and  v.  If  we  retain  only  terms  up  to  first 
order  (as  we  will  do  throughout  this  paper)  we  will  have 


n+1 


e(u 


i  V  ) 

n 


4 

u  ■  u 


(u 


n+1 


-  u  ) 

n 


'v 


n+1 


Vn) 


n+1 


V  .  ■  h(u  ,V  )  + 

n+1  ir  n  3u 


(u 


n+1 


\  'Hi 

-  U  )  +  r- 

n  ^v 


(v 


n+1 


v  ) 
n 


2S  +2K, 

u  2  u  3  n+1 
n  n 


(3) 


It  is  important  to  notice  that  Eq  (Q)  are  linear  in  u  and  v  . 

n+1  n+1 

We  will  make  use  of  the  property  of  linear  equations  that  the  solution 
is  the  sum  of  the  particular  solution  plus  the  sum  of  weighted  homo¬ 
geneous  solutions. 

If  ve  knew  uQ(t)  and  vn(t)  we  would  be  in  a  position  to  find 
un+^(t)  and  vnyj_(t)-  We  start,  then,  by  assuming  a  solution  which 
is  compatible  with  the  boundary  conditions  given  in  Eq  (7) 


uQ(t)  -  c 
v0(t)  -  o 


(9) 


Using  these  functions  to  evaluate  uq( t )  and  vQ(t)  at  every 
point  we  will  proceed  to  finding  u^(t)  and  v^(t)--our  new  estimate 
of  the  solution.  Begin  by  finding  the  particular  solution  corresponding 


to  the  initial  conditions 


u1(o)  *  o 

(10) 

v^(o)  a  o 


Solving  Eq  (8)  subject  to  these  conditions  yields 


U^(t)  =  qjt) 

•^(t)  -  q2U) 


(11) 


We  next  find  the  homogeneous  solutions  by  letting  all  terns  involving 
only  uq  and  vq  be  zero.  Solve  Eq  (8)  tvice  by  first  letting 


ux(o)  =  1 
v1(o)  *  o 


(12) 


and  second  letting 


u  (o)  a  o 
1 

v^o)  *  1 


(13) 


The  initial  condi tons  in  Eq  (12)  give  rise  to  the  solutions 


uxnl(t)  =  rx(t) 
v^Kt)  -  r2(t) 


(14) 


While  the  second  set  of  initial  conditions,  Eq  (13),  has  as  its  solution 
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u/l2(t)  -  S^t) 
vxh2(t)  -  e2(t) 

The  total  solution  Is  then 

Mt)  »  q1(t)  +  c1r1(t)  +  CgB^t) 
vl(t)  "  q2(t)  +  cxr2( t)  +  cgB2(t) 


(15) 


(16) 


where  the  weighting  coefficients  are  determined  by  using  the  boundary 
conditions  of  Eq  (9)-  We  may  now  continue  In  a  similiar  fashion  to 
find  u„(t)  and  v2(t),  and  u  (t)  and  v^(t),  etc.  until  we  are  satisfied 
that  further  computation  is  unnecessary.  Since  Quasilinearization  has 
the  property  of  quadratic  convergence- -each  iteration  has  the  effect 
of  virtually  doubling  the  number  of  correct  digits--very  few  iterations 
are  required. 

There  are  two  disadvantages  to  this  method.  The  first  is  that 

solving  for  the  weighting  coefficients,  involving  as  it  does  the  solution 

of  linear  equations,  can  lead  to  errors  if  the  matrix  of  the  homogeneous 

2 

solutions,  evaluated  at  the  requisite  points,  Is  ill-conditioned. 

The  second  is  that  the  solution  to  the  n+lst  solution  requires  that 

the  entire  nth  solution  be  stored  in  the  memory  of  the  computer.  This 

3 

can  quickly  lead  to  an  unacceptable  condition.  A  way  to  avoid  this 
is  to  store  only  the  initial  conditions  of  the  1st,  2nd,  ...,  nth 
solution  and  to  solve,  at  the  n+lst  iteration,  all  of  the  preceeding 
n  solutions  simultaneously.  We  still  have  a  positive  gain  for  we  have 
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t  railed  Limited  computer  storage  capability  for  increased  computation- - 
thc  thing  a  digital  computer  does  best. 

2.  DYNAMIC  PROGRAMMING 

We  begin  by  defining  an  o:timal--a  minimum  in  iur  ase--cost 
function 

f(c,T)  =  u°  JrV  =  *“un  F(u,u)dt  (17) 

Here  f(e,T)  represents  the  minimum  possible  cost  when  we  start  in  state 
c  with  tine  T  remaining  and  proceed  along  the  optimum  path  u(t).  It  is 
important  to  notice  that  if  we  have  no  time  (T  ■  0)  remaining  to  our 
process  we  can  do  nothing  that  will  change  the  value  of  the  incurred 
coot  regardless  of  the  state  of  the  system.  That  is  to  say 

f(c,o)  =  0  (18) 

The  principle  of  optimality  states,  'An  optimal  policy  has  the 
property  that  whatever  the  initial  state  and  initial  decision  arc,  the 
rci.iaining  decisions  must  constitute  an  optimal  policy  with  regard  to  the 
state  resulting  from  the  first  decision.11*1 

To  implement  the  principle  of  optimality  let  us  consider  the 
system  with  T  -  A  remaining.  The  cost  incurred  by  the  system  is,  to 
first  order,  F(c,v)\,  where  here  again  we  have 

v  -  u  (19) 

No  longer  are  we  in  state  c;  rather  ve  now  have 


u  ■  c  +  vA 


(20) 
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Applying  Eq  ( 17 )  vc  can  cay  that  the  minimum  (optimum)  coct 
a  procccn  initially  In  state  c  ♦  vh  with  T  -  L  remaining  ic 
T  -  ii).  If  wc  truly  wish  the  optimum  cost  function  ve  must 


Incurred  in 
f(c  +  v  , 
satisfy  the 


r elationchii 


f(c,T) 


rain 

v 


F(c,v)',  +  f(c  +  V 


(21) 


The  value  of  v  resulting  in  the  minimization  is  then  our  desired  v. 
Eq  (21)  can  now  be  exoanded  in  a  Taylor' s  series 


f(c,T)  -  “vN  F(  c,  v)i  +  f(c,T)  +  ~  vA  -  1 


(22) 


f(c,T)  appears  on  both  sides  and  is  Independent  of  v  so  it  can  be  can¬ 
celled.  We  can  substitute  in  for  F(c,v) 

F(c,v)  -  |v2  +  |  (23) 


Dividing  through  by  '  we  have 


0 


}f 


1 


(24) 


In  this  case  we  can  find  the  v  that  will  minimize  this  expression  analytically; 
generally  this  would  involve  search  techniques  on  a  computer.  Here  the 
optimum  v  is  found  using  basic  calculus  as 


v 


3f 

^c 


(25) 


Substituting  this  back  into  Eq  (24)  and  rearranging  leads  to 


3f  _  _  if  tt\2  +  K 
^T  2 \  Sc )  c 


(26) 
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Eq  (2b)  can  jc  solved  analytically  at  each  point  on  a  grid,  with  a 
suitable  number  of  joints,  in  the  c,T  plane.  At  each  point  we  would 
thou  iiave  botn  the  minimum  cost  to  be  incurred,  if  we  carried  out  the 
i  rooess  Li.  an  jptirnl  fashion  frot.  tiiat  point,  and  the  vaiuc  of  / 
corresponding  t.  the  optimal  jath  at  that  point.  I.,  other  words  we 
iiave  a  i'eeibn  >rocess  which  always  indicates  the  correct  direct!  n. 
t  rxec  i.-i.  Starting  at  any  point,  u(t)  is  then  o..jietcly  ..ocified 
n.  tens  .1  v(t). 

jj.  in  //d^in.r*  I!  wixDDIiIG 

Again  we  start  with  the  system  of  two  first  order  equations 
^  /on  ir,  Eq  (6)  subject  to  the  boundary  conditions  given  in  nq  ('/). 

■Ve  wish  to  imbed  oui’  original  problem  ii  a  more  general  one  valid  for 
any  initial  osition  d  and  any  initial  time  t.  It  is  obvious  that  the 
■orrect  choice  of  v  initially  will  be  a  function  of  initial  position 
a  .id  initial  tine,  Jc  note  this  dependence  by  saying; 

v(t)  *  r(d, t)  (27) 

II  w  we  -will  consider  the  process  at  tine  t  +  A. 
u( t+A)  ■  u(t)  +  v( t) \ 

*  d  +  r(d, t)A  ( 26) 

v(  t+A)  ■  r"d  +  r(d,  t)A,  t  +  A-1  “  v(T)  +  v  (t)A 

“  r(d,  t)  -  ^-A  (29) 


Expanding  the  left  side  of  Eq  (29)  we  have 
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r(d,T)  +  r(d,r)  -  r(d,T)  -  iU  (30) 

a 

Simplifying  and  rearranging  yields 

^  r^d,T^  ‘  ~2  (31) 

d 

Sq  (31)  represents  a  feedback  control  law  which  can  be  solved  for  any 
position  d  and  any  tine  r  subject  to  the  condition 

r(d,T)  •  0  (32) 

In  actually  implementing  Dynamic  Programming  and  Invariant 
Imbedding  digital  computers  would  be  used  to  numerically  solve  the 
feedback  control  laws.  The  accuracy  of  the  solution  then  becomes  a  function 
of  the  number  of  points  used  in  the  solution  grid  and  nay  represent  a 
computer  storage  problem.  Obviously  this  is  not  unique  to  the  methods 
discussed  but  rather  is  common  to  any  numerical  solution  of  partial 


differential  equations. 
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